Non-white noise and a multiple-rate Markovian closure theory for turbulence 
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Markovian models of turbulence can be derived from the renormalized statistical closure equa- 
tions of the direct- interaction approximation (DIA). Various simplifications are often introduced, 
including an assumption that the two-time correlation function is proportional to the renormalized 
infinitesimal propagator (Green's function), i.e. the decorrelation rate for fluctuations is equal to the 
decay rate for perturbations. While this is a rigorous result of the fluctuation-dissipation theorem 
for thermal equilibrium, it does not necessarily apply to all types of turbulence. Building on previous 
work on realizable Markovian closures, we explore a way to allow the decorrelation and decay rates 
to differ (which in some cases affords a more accurate treatment of effects such as non- white noise) , 
while retaining the computational advantages of a Markovian approximation. Some Markovian 
approximations differ only in the initial transient phase, but the multiple-rate Markovian closure 
(MRMC) presented here could modify the steady-state spectra as well. Markovian models can be 
used directly in studying turbulence in a wide range of physical problems (including zonal flows, 
of recent interest in plasma physics), or they may be a useful starting point for deriving subgrid 
turbulence models for computer simulations. 

PACS: 47.27.Eq, 47.27.Sd, 05.40.-a 
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I. INTRODUCTION 

Our derivation builds on and closely feJlows the work 
by Bowman, Krommes, and OttavianiEJ (we will fre- 
quently refer to this paper as BKO), on realizable 
Markovian closures derived from Kraichnan's direct- 
interaction-approximation (the DIA). The DIA is based 
on a renormalized perturbation theory and gives an 
integro-differential set of equations to determine the two- 
time correlation function. The DIA involves time inte- 
grals over the past history of the system, which can be 
computationally expensive. Markovian approximations 
give a simpler set of differential equations that involve 
only information from the present time. They approx- 
imate two-time information in the correlation function 
and in the renormalized Green's function by a decorrela- 
tion rate parameter. The structure of the equations we 
derive here is similar to the realizable Markovian closure 
(RMC) of BKO,lil but with extensions such as replacing a 
single decorrelation rate parameter with several different 
nonlinear rate parameters, to allow for a more accurate 
model of effects such as non-white noise. (As will be 
discussed more below, the RMC does include more non- 
white-noise effects than one might think at first.) 

The basic issue studied in the present paper can be 
illustrated by a simple Langevin equation (which will be 
discussed in more detail in the next section): 

(|+'?)vXi) = /W, (1) 

where rj is the decay rate and / is a random forcing or stir- 
ring term (also known as noise). As is well known, if / is 
white noise, then the decorrelation rate for tp is given by 
ry, so that in a statistical steady state the two-time corre- 
lation function {ip{t)ilj*{t')) = Cq exp(— 77|i — i'|) (assum- 
ing constant real rj here). However, if f(t) varies slowly 
compared to the I/77 time scale, then the solution to the 
Langevin equation is just ■ip{t) « f{t)/ r], and the decorre- 
lation rate for ijj is instead given by the decorrelation rate 
for /. Note that the Green's function (the response to a 
perturbation at time t') is still exp{—'ri{t — t')). Previous 
Markovian closures employed some variant of an ansatz, 
based on the fluctuation-dissipation theorem, that the 
two-time correlation function and the Green's function 
were proportional to each other. This is a rigorous re- 
sult for a system in thermal equilibrium, but may not 
necessarily apply to a turbulent system. The purpose of 
the present paper is to explore an extended Markovian 
closure, which we will call the Multiple-Rate Markovian 
Closure (MRMC), that allows the decorrelation rate of 
"0 to differ from the decay rate rj. 

In practice, the corrections due to non-white-noise ef- 
fects may be quantitatively modest, as the decorrelation 
rate for the turbulent noise / that is driving tp at a par- 
ticular wave number k is often comparable to or greater 
than the nonlinear damping rate 77 at that k. This is 
because the turbulent noise driving mode k arises from 



the nonlinear beating of other modes p and q such that 
k = p + q. Thus \p\ or |q| has to be comparable to 
or larger than |fc|, and will thus have comparable or 
larger decorrelation rates, since the decay rate rj is usu- 
ally an increasing function of |fe|. Furthermore, there 
are some offsetting effects due to the time-history inte- 
grals in the DIA's generalized Langevin equation that 
might further reduce the difference between the decay 
rate and the decorrelation rate. Indeed, past compar- 
isons of the RMC with the full DIA or with the full non- 
linear dynamics havei-geperally found fairly good agree- 
ment in many cases, tJa includiiur two-field Hasegawa- 
Wakatani drift-wgve turbulenceErQ and galactic dynamo 
MHD turbulence.Q Some of the results in this paper help 
to give a deeper insight into why this agreement is of- 
ten fairly good, despite the arguments of the previous 
paragraph, i.e., why the fluctuation-dissipation ansatz 
is often a reasonable approximation even out of ther- 
mal equilibrium. But there may be some regimes where 
the differences are important and the improvements sug- 
gested here would be welcome. These might include in- 
clude plasma cases where the wave dynamics can make 
rj vary strongly with the direction of k in some cases 
(with strong Landau damping in some directions and 
strong instabilities in other directions, for example), or 
non-steady-state cases involving zonal flows exhibiting 
predator-prey dynamics. 

Markovian closures such as the test-field model 
(TFM) or Orszag's eddy-damped quasinormal Markovian 
(EDQNM) closure have been extensively used to study 
turbulence in incomuressible fluids and plasmas. The 
introduction of BKOEl provides useful discussions of the 
background of the DIA and Markovian closures, and we 
will add just a few remarks here (there are also many 
reviews on these topics, such as Refs. |,|-@. The RMC 
developed in BKQEI is simflar to the EDQNM, but has 
features that ensure realizability even in the presence of 
the linear wave phenomena exhibited by plasmas (e.g. 
drift waves) and rotating planetary flows (e.g. Rossby 
waves). "Realizability" is a property of a statistical clo- 
sure approximation that ensures that, even though it is 
only an approximate solution of the original equations, 
it is an exact solution to some other underlying stochas- 
tic equation, such as a Langevin equation. The absence 
of realizability can cause serious physical and numerical 
problems, such as the prediction of negative or even di- 
vergent energies. The RMC reduces to the DIA-based 
version of the EDQNM in a statistical steady state, so 
in some cases the issue of realizability is only important 
in the transient phase as a steady state is approached or 
in freely decaying turbulence. Realizability may also be 
important in certain cases of recent interest among fusion 
researchers where oscillations may occur between various 
parts of the spectrum (such as predator-p rev t ype oscil- 
lations between drift waves and zonal flowaljllil) where a 
simple statistical steady state might not exist, or where 
one is interested in the transient dynamics. Unlike some 
Markovian models that differ only in the transient dy- 
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namics, the Multiple-Rate Markovian Closure presented 
here could also alter the steady-state spectrum. 

Our results apply to a Markovian approximation of 
the DIA for a generic one-field system with a quadratic 
nonlinearity. They are immediately applicable to some 
simple drift-wave plasma turbulence problems, Rossby- 
wave problems, or two-dimensional hydrodynamics. Fu- 
ture work could extend this approach to multiple fields, 
similar to the covariant multifield RMC of BKOEl or their 
later realizable test-field model.H Multiple field equations 
can get computationally difficult (with the compute time 
scaling as n^, where n is the nurnher of fields), though 
two-field studies-, have been doneO and disparate scale 
approximationala or other approximationaiS might make 
them more tractable. In addition to their direct use 
in studying turbulence in a wide range of systems, the 
Markovian closures discussed here might also be use- 
ful in deriving subgrid turbulence models for computer 
simulations .EZnlfl 

While our formulation is general and potentially ap- 
plicable to a wide range of nonlinear problems involving 
Markovian approximations, we were motivated by some 
recent problems of interest in plasn3|a_p|liysics and fusion 
energy research, such as zonal flows.Ej^EHl Initial analytic 
work elucidating the essentials of nonlinear zonal, flow 
generation used weak-turbulepee approximationsESEJ or 
secondary-instability pialysisEj Recent interesting work 
by Krommes and KimEj uses a Markovian statistical the- 
ory to extend the study of zonal fiows to the strong tur- 
bulence regime. An important question is why the strong 
generation of zonal flows seen near marginal stability is 
not as important in stronger instabilitji,r|egimes (i.e., why 
is the Dimits nonlinear shift finite?). ca~E3 A strong tur- 
bulence t heory is needed to study this. An alternative 
approach ,EllE3 which has been fruitful in providing the 
main answers to the finite Dimits shift question, is to 
analyze the secondary and tertiary instabilities involved 
in the generation and breakup of zonal flows. That work 
suggests that a complete strong-turbulence Markovian 
model of this problem would also need multi-fleld and 
geometrical effects (involving at least the potential and 
temperature fields, along with certain neoclassical effects 
in toroidal magnetic field geometry). 

Based on the reasoning immediately following Eq. (|l|) 
above, one might think that the assumption that the two- 
time correlation function and the Green's function are 
proportional to each other is rigorous only in the limit 
of white noise (which has an infinite decorrelation rate). 
The Realizable Markovian Closure has been shown to 
correspond exactly to a simple Langevin equation (where 
the effects of the turbulence appear in nonlinear damping 
and nonlinear noise terms), for which this might appear 
to be the implication. However, the mapping from statis- 
tically averaged equations (such as Markovian closures) 
back to a stochastic equation for which it is the solution, 
is not necessarily unique. In particular, the full DIA cor- 
responds to a generalized Langevin equation (Eq. ( ^ ) 
below), in which the damping term r]ip{t) in the sim- 



ple Langevin equation is replaced by a time-history in- 
tegral operator. As we will find, it is then possible for 
the two-time correlation function and the Green's func- 
tion to be proportional to each other even when the noise 
has a finite correlation time. This allows the fiuctuation- 
dissipation theorem (which is rigorous in thermal equi- 
librium) to be satisfied without requiring the noise to be 
white (since the noise is not necessarily white in thermal 
equilibrium). Thus the fluctuation-dissipation ansatz of 
BKO is a less restrictive assumption than one might have 
at first thought. [It should be noted that previous Marko- 
vian models account implicitly for at least some non- 
white-noise effects. For example, in the calculation of 
the triad interaction time 9kpq = i/ivk + Vp + Vg) fo^' 
three- wave interactions, finite values of the assumed noise 
decorrelation rate rjp + r]q are used.] 

Nevertheless, there is still no reason in a general sit- 
uation that the two-time correlation function and the 
Green's function be constrained to be proportional to 
each other. As described elsewhere, there may be regimes 
where the resulting differences between the decorrelation 
rate and the decay rate are significant. 

The outline of this paper is as follows. Sec. (0) 
presents some of the essential ideas of this paper for a 
very simple Langevin equation, and includes a section 
motivating the choice of the limit operator introduced by 
BKCB to ensure realizability. Sec. (MW) presents a more 
detailed calculation of the non-white Markovian model 
for a simple Langevin equation, including the effects of 
complex damping rates (to represent the wave frequency) 
and the issue of Galilean invariance. The model is com- 
pared with exact results in the steady-state limit, and 
then an extension to time-dependent Langevin statistics 
is presented (along with, in Appendix |a| a n alternative 
proof of realizability for this case). Sec. (IV) presents the 
notation of the full many-mode nonlinear equations we 
will solve and summarizes the direct interaction approxi- 
mation (DIA), which is our starting point. Sec. (M) sum- 
marizes how the non-white Markovian approach is de- 
rived for the steady-state limit (with further details given 
in Appendix ^), while Sec. (VI) presents the full non- 
white Markovian approximation for the time-dependent 
DIA. Sec. (VII) discusses some important properties of 
these equations, including the limits of thermal equilib- 
rium and inertial range scaling, and some difficulties due 
to the lack of random Galilean invariance in the DIA and 
described in Appendix ^ The conclusions include some 
suggestions for future research. 



II. SIMPLE EXAMPLES BASED ON THE 
LANGEVIN EQUATION 

Here we expand upon the analogy given in Sec. (|) us- 
ing a simple Langevin equation, which provides a useful 
paradigm for understanding the essential ideas consid- 
ered in this paper. Since realizable Markovian closure 
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approximations to the DIA can be shown to correspond 
exactly to an underlying set of coupled Langevin equa- 
tions, the analogy is quite relevant. In this section we 
will consider heuristic arguments based on some simple 
scalings; later sections will be more rigorous. 
Consider the simple Langevin equation 



d_ 

dt 



(2) 



where 77 is a damping rate and /* is a random forcing 
or stirring term (also known as "noise"). [Here we now 
force Eq. (|l|) with the complex conjugate of /, for con- 
sistency with the form of the equations used later for a 
generic quadratically nonlinear equation.] The statistics 
of the noise are given by a specified two-time correla- 
tion function Cf{t,t') = {f{t)f*{t')). [In the whitc-noisc 
limit, Cf{t,t') = 2DS{t — t'), and the power spectrum of 
the Fourier-transform of f{t) is independent of frequency, 
and is thus called a "white" spectrum.] The Langevin 
equation is used to model many kinds of systems exhibit- 
ing random- walk or Brownian motion features. Here we 
can think of ip as the complex amplitude of one compo- 
nent of the turbulence with a specified Fourier wave num- 
ber k. Note that rj may be complex (representing both 
damping and wave-like motions) and represents both lin- 
ear and nonlinear (renormalized) damping or frequency 
shifts due to interactions with other modes. The random 
forcing /* represents nonlinear driving by other modes 
beating together to drive this mode. 

The response function (or Green's function or propa- 
gator) for this equation satisfies 



d 



(3) 



which easily yields R{t, t') = exp(-ry(i - t'))H{t - t') if 
r] is independent of time, where H{t) is the Heaviside 
step function. The solution to the Langevin equation 
is just %l){t) = f^dtR{t,t)f*{t) (for the initial condition 
■0(0) = 0). It is then straightforward to demonstrate the 
standard result that, if / is white noise and the long- 
time statistical steady-state limit is considered, then the 
correlation function for is 

C{t,t') = m)r{t')) = Coexp(-rKt- i')) 

for t > t' , where Cq — 2D/{r] + ij*) (we emphasize def- 
initions with the notation =). [For t < t' , one can use 
the symmetry condition C{t,t') — C*{t',t).] That is, 
the decorrelation rate for -0 is just 77. This is equivalent 
to the assumption in a broad class of Markovian models 
that the decorrelation rate for ip is the same as the decay 
rate of the response function. 

However, consider the opposite of the white-noise limit, 
where / varies slowly in time compared to the l/rj time 
scale. Then the solution of Eq. (||) is approximately 
ip{t) — f*{t)/ri, and the decorrelation rate for ^p will 
be the same as the decorrelation rate for /*. In this 



limit, the assumption in many Markovian models that 
the decorrelation rate is rj is not valid. 

Denote the decorrelation rate for /* as r/j, and the 
decorrelation rate for ip as rjc- Then one might guess 
that a simple Pade-type formula that roughly interpo- 
lates between the white- noise limit 77/ S> and the op- 
posite "red-noise" limit r]f <^rj would be something like 
l/ric ~ l/?7-f 1/77;, or 



(4) 



In fact, we will discover in the next section that more 
detailed calculations give similar results in the limit of 
real 77 and 77/, though the formulas are more complicated 
in the presence of wave behavior with complex 77 and 77/. 

We note that in many cases of interest, the noise decor- 
relation rate rjf turns out to be of comparable magnitude 
to 77 (for example, if the dominant interactions involve 
modes of comparable scale). In this case, while the white- 
noise approximation is not rigorously valid, the correc- 
tions to the decorrelation rate considered in this paper 
might turn out to be quantitatively modest, ~ 50%. Fur- 
thermore, in the case of the full DIA and its correspond- 
ing generalized Langevin equation, we will find additional 
corrections that can, in some cases, offset the effects in 
Eq. (Q) and cause -qc to be closer to 77. 

Before going on to the more detailed results in the next 
section, we considjer the meaning of an operator intro- 
duced in the BKOEl derivations in order to preserve real- 
izability in the time-dependent case, where rj{t) varies in 
time and may be negative (transiently), representing an 
instability. [In order for a meaningful long-time steady- 
state limit to exist, the net 77 (which is the sum of lin- 
ear and nonlinear terms) must eventually go positive to 
provide a sink for the noise term. But it is important 
to preserve realizability during the transient times when 
77 may be negative.] Based on arguments about symme- 
try and the steady-state fluctuation-dissipation theorem, 
they initially proposed a time-dependent ansatz of the 
form 



exp 



dt r]{t) 



(5) 



(for t > t'), where C(t) = C{t,t) is the equal-time 
covariance. Later in their derivation, they state that 
in order to ensure realizability, 77(f) in this expression 
had to be replaced with V{rj{t)), where the operator 
7^(77) = Jier}H{Kerj) + 7lm77 prevents the real part of 
the effective 77 in Eq. (H) from going negative. 

Physically this makes sense for the following reasons. 
Consider Eq. (^) with white noise / (thus ignoring the 
non- white-noise effects). Then in a normal statisti- 
cal steady state where r]{t) is constant and Re 77 > 0, 
Eq- (^ properly reproduces the usual result C{t,t') — 
Co exp(— 77|i — t'\). However, if Re 77 < (which it might 
do at least transiently in the full turbulent system con- 
sidered later), then Eq. (0) can't reach a steady state. 
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and the solution is eventually ip{t) = -ipo exp(— ryt) = 
exp(— i Imryi), after an initial transient phase. 

Thus Cit,t') = Ci/2(t)Ci/2(t')exp(-iIni?7(t - t')), in 
agreement with and providing an additional intuitive ar- 
gument for BKO's modified form of Eq. (||), including 
the V{ri) operator. [There may be an initial phase where 
the noise term / in Eq. (||) dominates and causes C(t) 
to grow Hnearly in time, C(t) = {tp{t)%p* {t)) = 2Dt. But 
eventually the unstable "qip term will become large enough 
to dominate and lead to exponential growth of ip] 

The model we will introduce below replaces 77 in Eq. (j^) 
with a separate parameter 77c, and develops a formula to 
relate 77c to other parameters in the problem such as 
77 and 77/. In the white-noise limit, the formula for rjc 
automatically reproduces the effects of the V limiting 
operator, as will be described in the next section and 
in Appendix (|^. But numerical investigation of non- 
white noise with wave dynamics (Im77 7^ or Im77/ 7^ 0) 
uncovered cases where the V limiting operator is still 
needed to ensur e real izability. This will be explained at 
the end of Sec. ( UlICD . 

We considered naming the method described in this 
paper the Non- White Markovian Closure since, for the 
simple Langevin equation considered here and in the next 
section, the decorrelation rate and the decay rate are 
equal only in the white- noise limit, and this approach al- 
lows these rates to differ. [Alternatively, to emphasize 
the flexibility of this method one might have called it the 
Colored-Noise Markovian Closure since instead of being 
restricted to white-noise (a uniform spectrum) we can al- 
low a noise spectrum of width Slo Her/f peaked near 
an arbitrary frequency uj Imrif. In other words, this 
closure can model spectra with a range of possible col- 
ors.] However, as we will discuss further, while a simple 
Langevin equation is sometimes used to demonstrate re- 
alizability of Markovian approximations, the DIA is actu- 
ally based on a generalized Langevin equation involving 
a non-local time-history integral (compare Eq. (^) with 
Eq. (^o|)). Because non- white fluctuations enter not only 
by making the noise term non-white but also by affecting 
this time-history integral, it is possible for the decay-rate 
and the decorrelation rate to be equal even in some cases 
where the noise is not white (as indeed is the case in ther- 
mal equilibrium where the fluctuation-dissipation theo- 
rem must hold but the noise is not necessarily white). We 
thus favor the name Multiple-Rate Markovian Closure 
(MRMC), to emphasize that the method developed here 
is a generalization of the previous Realizable Markovian 
Closure (RMC) to allow for multiple rates (i.e., separate 
decay and decorrelation rates). 



III. DETAILED DEMONSTRATION OF THE 
MULTIPLE-RATE MARKOVIAN METHOD 
WITH THE LANGEVIN EQUATION 



equation. The steps in the derivation are quite similar 
to the steps that will be taken in the following sections 
for the case of the more complete DIA for more compli- 
cated nonlinear problems, and thus help build insight and 
familiarity. In this section, we will be introducing vari- 
ous approximations that may seem unnecessary for the 
simple Langevin problem, which can be solved exactly 
in many cases (for simple forms of the noise correlation 
function). But these are the same approximations that 
will be used later in deriving Markovian approximations 
to the DIA, and so it is useful to be able to test their 
accuracy in the Langevin case. 

Our starting point is the Langevin Eq. (^, but we 
allow ri(t) to be a function of time, so that the solution 
to Eq. (0) for the response function is 



i?(t, t') = exp 



dtf]{t) H{t-t') 



(6) 



(instead of the solution given immediately after Eq. (|^), 
which assumes that 77 is independent of time) . The solu- 
tion to the Langevin equation is 



dtR{t,t)f*{t). 



(7) 



In principle it is possible to calculate directly two-time 
statistics like C{t,t') = {tjj{t)ilj* {t')) from this, but in 
practice it is often convenient to consider instead the dif- 
ferential equation for dC{t,t')/dt, which from Eq. (|^) 
and Eq. (Q) is 



dtR*{t',t)C*j{t,t), 



(8) 



where the noise correlation function is defined as 
Cf{t,t') = {{f{t)f*{t')), and we have assumed that the 
initial condition ipiO) has a random phase. This equation 
is the analog of the DIA equations fo r th e two-time corre- 
lation function (compare with Eq. (39a) and Eqs. (^l|)), 
but with an integral only over the noise and no nonlinear 
modification of the damping term. 

We define the equal-time correlation function C{t) in 
terms of the two-time correlation function C{t,t') as 
C{t) = C{t,t) = {tp{t)ip*{t')) (note that these two func- 
tions are distinguished only by the number of arguments). 
Then 



dCjt) 
dt 



2Rcri Cit) = 2Rc J dt R*{t,t)C}{t,t). (9) 



This is the analog of the DIA equal-time covariance equa- 
tion, Eq. dH). 



In this section, we demonstrate the Multiple-Rate 
Markovian approach starting with a simple Langevin 
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A. Langevin statistics in the steady-state limit 

Consider the steady-state hmit where t, i' — > cx) (but 
with finite time separation t — t'), and assume the noise 
correlation function has the simple form Cf{t,t') — 
Cfoexp[—rif{t — t')] for t > t' . In this section we as- 
sume rj and r]f are time-independent constants. The 
response function reduces back to its steady-state form 
R{t,t') = exp[-7-i{t-t')]H{t-t'). Then Eq. (|) in steady 
state gives 



Co ^ lim C(t) - ^^/Oi,,(^)(^'|'^))Jj!^^p- 



(10) 



Writing rj = v+itu and rjf = Vf+icuf in terms of their real 
and imaginary components, and denoting the frequency 
mismatch Ato = uj + ujf (remember, because the complex 
conjugate /* is used as the forcing term, resonance occurs 
when Im(77) = Im(77p) this can be written as 



Co 



Cfo {t^ + vf) 
V (i/ + t^/)2 + (Aw)2 



(11) 



This has a familiar Lorentzian form characteristic of res- 
onances. 

To find the two-time correlation function, the time in- 
tegral in Eq. (||) can be evaluated for i > i' to give 



at J rj^ + rjj 



exp[~r,}{t~t')]. (12) 



With the steady-state boundary condition C{t = t' ,t') = 
Co, this can be solved to give 



C{t,t') = Co 
+ Co 



1 - 



exp[— 77(i — t') 



Re{ri + rif){r] - 7]*) 
Re{r]){r] + r]f) 

Re(, + ,,)(,-,;)^^P["^/^^-^)]- 



In the white-noise limit, \rif\ ^ \ri\, this reduces to the 
standard simple result C{t,t') = Co exp[-'i]{t — t')]. But 
in the more general case of non-white noise, the two- 
time correlation function is more complicated. [Despite 
the apparent singularity in the denominator, it is can- 
celled by the exponentials so that C{t, t') is well-behaved 
in the limit f] ^ V}-] In the context of the turbulent in- 
teraction of many modes, Cf{t,t') and thus C{t,t') may 
be very complicated functions. Even if the noise cor- 
relation function has a simple exponential dependence 
Cf{t, t') oc exp[—rif{t — t')], we see that the resulting cor- 
relation function for is more complicated. 

Consider the task of fitting this complicated C{t,t') 
with a simpler model of the form 



C,^od{t,t') = Coexph77c(i-t')] 



(14) 



(for t > t'). One way to define the effective decorrela- 
tion rate "qc might be based on the area under the time 
integral. 



dt C,„od(i, t ) = — 



dt'C{t,t'). (15) 



This can be evaluated either by directly substituting 
Eq. ([T^), or by taking a time average of Eq. (^2|); the 
same answer results either way. It turns out that in the 
later versions of this calculation it is easier to determine 
rjc by integrating the governing differential equation over 
time. Operating on Eq. (|l2|) with J*i^ dt' and using 



dt 



,dC{t,t') 



we find 



dt 
1 

lie 



d_ 

di 



dt'C{t,t') -C{t,t), (16) 



1 ^ Re(?7)(?? + ??/) 

7] Rc{7] + 7]f)7]7]j 



(17) 



This recovers the white-noise limit rjf ^ rj and the red- 
noise limit rjf <^ rj discussed in Sec. ||. In the limit of 
real rj and real 77/ it simplifies to the Pade approxima- 
tion 7^c = VV f / iv + V f) 8.1so suggested in the introduc- 
tion. However, there is a problem with Eq. ( p7| ) related to 
Galilean invariance. Suppose we make the substitutions 
Tp = 'ibeKp[iLU2t] and /* = /* exp[iLi;2t] into the Langevin 
Eq. (H). Then it can be written as 



|+r})^ = r(0. 



(18) 



where fj = rj + 1^2, and the results should be the 
same if written in terms of the transformed variables. 
In particular, the correlation function should trans- 
form as {7p{t)7P*{t')) = cxp[-iuj2{t - t')]{7p{t)'tl;*{t')) = 
eyip[—iuj2it ~ t')]C{t, t'). Thus the decorrelation rate f}c 
for Tp should be related to the decorrelation rate 77c for 7p 
by fic = Vc + i^i- The decorrelation rate for the trans- 
formed noise term /* also transforms as 77^ = 77J -I- iijJ2- 
In the case of fiuid or plasma turbulence where Tp rep- 
resents the amplitude of a Fourier mode oc exp[zfc • a;] 
and /* represents the amplitude of two modes with wave 
numbers p and q beating together to drive the fc mode 
(so p + q = k), these transformations correspond to a 
Galilean transformation to a moving frame x = Xo + vt, 
with UJ2 — k ■ V. 

So all results should be independent of lu2 under the 
transformation 7] = 7) — 10^2, rj} — ~ *^2, (thus 
Vf — Vf + ^^2), Vc = Vc ~ i^2- Eq- (|ll| ) satisfies this, 
but Eq. (^^ fails this test. This problem and its-jSolu- 
tion is described in the review paper by KrommesJHj who 
shows it is related to other differences in various previous 
Markovian closures. The problem can be traced to the 
definition of Eq. (|l5|), which doesn't satisfy the invari- 
ance for general forms of C(i, t'). For example, we could 
have multiplied the integrand in Eq. (Ha) by an arbitrary 
weight function (such as exp[— zu;2(t — 7)]) before taking 
the time average, and the results would have changed. 
The way to fix this problem is to do the time average in 
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a natural frame of reference for ^ that accounts for its 
frequency dependence. This leads us to the definition: 



'-'0 



dt' Cl,^{t,t')C{tX) 



(19) 



This corresponds to fitting Cmod{t,t') to C{t,t') by re- 
quiring that both effectively have the sarpe projection 
onto the function Cmodit, t'). [As Krommescj points out, 
using the invariant definition Eq. (^9|) instead of Eq. ( p5| ) 
is a non-trivial point needed to ensure realizability and 
avoid spurious nonphysical solutions in some cases.] 

Operating on Eq. ( p^) with dt' C^^^{t, t'), using a 
generalization of Eq. (tt6|), and doing a little rearranging 
yields 



Vc 



Cfoir]c + Vc) 



(20) 



This is properly invariant to the transformation described 
in the previous paragraph. Solving for rjc while leaving 
77^ on the other side of the equation, eventually leads to 



Vc = 



vv*f Rc(?? + Vf) + iv*c ^^{yy}) 

(77 + Tj*c) Re(77 + Tjf) + (77; + Tj*) Re(77/) ' 



(21) 



If we consider the limit where ??, ?7/, and thus rjc are all 
real, this simplifies to the form 



Vc 



7777/ 



V + Vf +VC 



(22) 



This is similar to (but more accurate than) the rough in- 
terpolation formula Eq. (^ suggested in the introduction. 
This kind of recursive definition, with rjc appearing on 
both sides, is a common feature of the steady-state limit 
of theories based on the renormalized DIA equations, and 
can be solved in practice by iteration, or by considering 
the time-dependent versions of the theories. In Eq. ( |2^ ) 
with real coefficients, one can easily solve this equation 
for rjc , but the solution is much more difficult in the case 
of complex coefficients in Eq. (^l|) . The resulting calcula- 
tion is laborious, so we used the symbolic algebra pack- 
age MapleKl to solve for rjc with complex coefficients. 
Looking at the real and imaginary parts of Eq. ( ^l|) sep- 
arately eventually leads to a quadratic equation and a lin- 
ear equation to determine the real and imaginary parts 
of i]c- Unfortunately it takes 16 lines of code to write 
down the resulting closed-form solution (though perhaps 
there are common subexpressions that would simplify it). 
(Maple worksheets that show this calculation and check 
other main results in this paper are available online.Ell) 
This is tedious for humans but easy to evaluate in For- 
tran, C, or other computer language. On the other hand, 
this is only helpful for the simple Langevin problem any- 
way since direct solution is not really practical for the 
full nonlinear problem considered by the DIA, where the 
noise term of the Langevin equation is replaced by a sum 



over many modes. In many cases of interest, the noise 
decorrelation rate 77/ turns out to be comparable in mag- 
nitude to r/, so iteration of Eq. ( ^l| ) usually converges 
quickly. (However, there are limits where convergence 
is slow, such as some strongly non-resonant cases where 
Re?7 is very close to Re 77/ and both are very small com- 
pared to lm(r]f + 77).) The other option is to consider 
the time-dependent problem, the topic of the subsection 
after next, which effectively performs an iteration in time 
as a steady state is approached. 



Re C(t,t') 




exact 
model 
wn 



Ti,= 0.25 



FIG. 1. KeC{t,t')/Co vs. t — t', for the exact Langevin re- 
sult of Eq. ( p^ ) , for the Multiple-Rate model with decorrela- 
tion rate 77c given by Eq. (^l|) , and for the simple white-noise 
assumption C{t,t') — Coexp(— 77|t — t'\). Time is normalized 
such that 77 = 1, and the value of 77/ is noted in each figure. 
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exact 
model 
wn 



FIG. 2. Real and imaginary parts 

for the same three functions as 
rif = 0.25 - 4i. Note that ImC = 
in this and later figures. 

1 



Tif= 0.25 - 4 i 



oiC{t,t')/Co vs. t-t', 
in Fig. 1, but with 
for the white-noise case 




exact 
model 
wn 




FIG. 4. Real and imaginary parts of C{t,t')/Co vs. t — t', 
for the same three functions as in Fig. 1, but with r]f = l — 4i. 



Re C(t,t') 



exact 
model 
wn 




FIG. 3. Real and imaginary parts 
for the same three functions as in FIe 



oiC{t,t')/Co vs. t-t', 
;. 1, but with rif = 1 — 



FIG. 5. Real and imaginary parts of C{t,t')/Co vs. t — t', 
for the same three functions as in Fig. 1, but with rjf = 1— 16i. 
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Re C(t,t 



Im C(t,t' 




-0.4 
-0.6 
-0.8 



-1J 



■ exact 
' model 
wn 



11,= 4- 16 i 



FIG. 6. Real and imaginary parts of C{t,t') /Co vs. t — t' , 
for the same three functions as in Fig. 1, but with r/f = 4— 16i. 



arise when Imrjf ^ 0, while the multiple-rate model 
does a fairly good job of capturing the real and imag- 
inary parts of C{t,t') in most cases. The most challeng- 
ing case for even the multiple-rate model is depicted in 
Fig. (5), where there is a large frequency mismatch but 
comparable decorrelation rates, Re?7/ ~ Rery. However, 
Eq. (|ll|) shows that the amphtude, Cq ~ 2C/o/(Aw)^ '-^ 
2C/o/(Ini(?7 — 'n*f)Y 1 '^ill be small in this strongly non- 
resonant case, and perhaps does not matter much com- 
pared to resonant interactions in realistic many-mode 
turbulence cases. Strongly non-resonant cases are eas- 
ier to model with disparate values of Re 77/ and Re 77, as 
shown in Fig. (6) and Fig. (2), because interference effects 
are less important. To do better for the non-resonant case 
with Re?7/ ~ Rery would probably require a more elab- 
orate two-exponential model than Eq. ([l^), to allow for 
the constructive and destructive interference effects rep- 
resented in Fig. (5). Of course, for the simple Langevin 
case of this section, such a model could exactly repro- 
duce Eq. (p^, although for more complicated cases it 
would again become a model to be fit to the true C{t, t') 
dynamics. (Another approach, which might improve the 
long-time fit a bit, might be to use C^^^{t, t'){t—t') as the 
weight function in Eq. ( p^ ) instead of just C'^^^^{t,t').) 



C. Time-dependent Langevin statistics 



B. Comparison of the Multiple- Rate model with 
exact Langevin result 

Figs. (1-6) provide a comparison of the exact and 
model results for various parameters. The exact 
Langevin solution for C{t,t')/Co is given by Eq. (13). 
The curves labeled "model" are for the Multiple-Rate 
Markovian model Cmod{t,t')/Co — exp{—r]c\t — t'\), 
where rjc is obtained by solving Eq. (pi|). The curves 
labeled "wn" are the results for a simple white-noise as- 
sumption C{t,t')/Co = exp{—ri\t — t'\). The plots show 
both the real and imaginary parts of C{t, t'), except when 
C{t,' t) is purely real. 

The results are shown in Figs. (1-6) for a variety 
of parameters. We choose 77 = 1 as a standard nor- 
malization in all cases. Only the frequency mismatch 
{Alo = lm(ri — ryp) between the oscillator and the ran- 
dom driving term and the relative decorrelation rate 
(Re(?7y)/ Re(?7)) can matter. Thus we choose a frame 
of reference where Im 77 = and any frequency mismatch 
is reflected in the value of the noise frequency Im77/. 

These comparisons show that the non-white-noise 
Multiple-Rate model for rjc does fairly well in most cases. 
All formulas of course agree well in the white-noise limit 
of Re 77/ 3> Re ij. The errors of the white- noise model are 
particularly large in the "red-noise limit" Re 77/ <C Re 77, 
though they are noticeable even if Re 77/ ^ Re 77. The 
white-noise model has a purely real correlation func- 
tion in all cases, thus missing the frequency shifts that 



We now return our attention to the more general 
Langevin problem with time-dependent 77 (t) and time- 
varying statistics for the noise term f*{t). That is, for 
generality, we also allow the noise amplitude (given by 
the equal-time covariance Cf{t) = Cf{t, t)) and the noise 
decorrelation rate to vary in time. Our choice of a self- 
consistent model for Cf{t,t') to accomplish this is moti- 
vated by BKO's demonstration that the following form 
is a realizable correlation function: 



Cf{t,t) = Cy^(i)exp 



Cy\t) (23) 



(for t > t). [BKO show this is realizable as long as 
Re{r]f{t)) > almost everywhere.] Using this expres- 
sion, Eq. (|^) can be written as 



dCjt) 
dt 



2Ref]{t)C{t) = 2Re Cy^{t)e*{t), (24) 



where 



e{t)= J dtRit,t)exp dt"rif{t") C]''^{1). (25) 

Taking the time derivative of this expression, and using 
Eq. (^), leads to 



dQ{t) 
dt 



-m+mit)mt)+c]'\t), (26) 



which is more convenient to use in a time-dependent cal- 
culation than Eq. (p5|). The initial condition is 0(0) = 0. 
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Eq. ( |24D and Eq. (g6|) can be used to determine the 
equal-time covariance C{t), but how can we determine 
the decorrelation rate rjc from the two-time correlation 
fimction C{t,t')? [In the full nonlinear equations used 
for the DIA, for one mode appears in noise terms for 
other modes, and so we would like to know the decorre- 
lation rate as well as the amplitude C{t).] Even in the 
steady-state limit of the previous section, we found that 
the full two-time correlation function C{t, t') had a more 
complicated form than a simple exponential, and so we 
fit a simpler model C^o^{t, t') to it in order to determine 
an effective decorrelation rate r]c- 

We follow a similar procedure here. We again use 
BKO's form for a realizable time-dependent two-time cor- 
relation function to provide a model of C{t,t'), 



cxp 



- dt"iic{t") 



(for t > t'). Consider the integral 



A{t) 



dt' C*^,S.t')C{t,t'). 



(27) 



(28) 



This is the time-dependent analog of Eq. (|19[). Rather 
than try to use this to determine "qc directly, it is more 
convenient to again take time derivatives. If C{t,t') in 
Eq. (^) is replaced with Cnwd{t,t') of Eq. (^, then 

^^(^) - c^t) - [^c{t) + + A. (29) 



dt 



C{t) dt 



If we instead calculate dA/dt with the full C{t,t') in 
Eq. (p8|), and use Eq. (||) to evaluate dC{t,t')/dt, then 



9Mt) ^2. N r . N . . M . 1 9C{t) ^ 



dt 
where 



2C{t) dt 



+ Ql{t)C'l'{t)Cy\t), 



dt' 



C^/^{t) 



r* _ c*f 

/ dtR*(t',t]^ 

Jo c] 



c}{t,i) 



(30) 



(31) 



Taking the time derivative of this, and using Eq. (|2^) for 
Cfit,t), gives 



^c'/^t)eit)~[^cit) + vfit)]e3it). (32) 



dt 

Equating Eq. ( ps| ) and Eq. (^0|), one can then solve 
for the effective decorrelation rate rjc ■ Using Eq. ( p^ ) to 
eliminate the dC(t)/dt term, the result is 



Vc{t) = V [v{t)-Rcrj{t) 



Cy^{t) ReQ{t) 



C{t) 



Qi{t)c^/^{t)cy\t) 

Mt) 



(33) 



where we have added the V operator to enforce realiz- 
ability for the reasons discussed below. Here V{z) = z 
if Rcz > and Viz) = ilmz if Rez < 0. Substituting 
Eq. (H) into Eq. (||) gives 



dAjt) 
dt 



C^t)-2Re{r^{t)+r^c{t))A{t) 

C'^^^(t) 
+ 2ReQ{t)^j^A{t). 



(34) 



Eqs. (|2^), (g6D, and (^2[]3^) provide a complete set of 
equations that can be integrated forward in time. They 
comprise a Markovian closure theory (including non- 
white noise effects) for the time-dependent Langevin 
equation. The relevant initial conditions are discussed 
below. This set of equations can be used to determine 
the amplitude C{t) and the effective decorrelation rate 
ijcit) used to model the two-time behavior C{t, t'). 

In a normal long-time statistical steady state, where 77, 
rff and C/ are constants (and Re(?7) > and Re(?7/) > 0), 
then one can show that the second and third terms on the 
right-hand side of Eq. ( |33| ) cancel and that it reproduces 
the steady-state result for rjc in Eq. (20). 

Consider the behavior of these equations in an unstable 
case, with Re 77 = —7 < 0. For simplicity, assume the 
coefficients r], rjf and Cf are all constant in time, with 
Re{r]f) > 0. Then one can show that C{t) eventually 
grows as exp(27f), while 8(t) ~ exp((7 — r]f)t) grows 
more slowly, so that the third term on the right-hand side 
of Eq. (|3^) vanishes. The fourth term on the right-hand 
side of Eq. ( ^ ) also vanishes because Q3 ~ exp(27t) 
while A ^ exp(47t). In this limit, tjc = f] — Re(77). 

Thus with constant coefhcients, the two cases of pos- 
itive or negative Re 77 will, at least in the long-time 
limit, naturally reproduce the limiting operator V(r]) =p. 
ReriH{Reri) + ilmrj, which was introduced by BKOEJ 
to preserve realizability for the assumed form of C{t,t') 
in Eq. (p7|). In the white- noise limit rjf 3> rj, it is 
straightforward to show that realizability is ensured for 
all time, not just in the long-time limit (see also Ap- 
pendix (^). These results might suggest that the V op- 
erator in Eq. (33) is not needed, if its argument always 
has a positive real part anyway. However, by numeri- 
cally integrating Eqs. (||), (|2|), and (|2[]3|), we have 
found cases where this is not true and the V operator is 
needed in Eq. ( |33| ) to enforce the realizability condition 
Re 77c > 0. [Without the V operator, Rct^c will tran- 
siently go negative in some strongly non-resonant cases 
such as 77 = 1 and 77/ — 0.25 -I- 16?.] Eqs. (24, ^6|) are an 



exact system of equations for the equal time covariance 
C{t) for Langevin dynamics, which ensures that C{t) 
always positive. But according to Theorem 2 of BKOl 
(and Appendix ^ of the present paper). Re rjc > is 
necessary for Cniod(t,t') as given by Eq. ( p7| ) to be a 
realizable two-time correlation function. This may be 
important if C'mod(t,t') is in turn used in a noise term 
driving some other Fourier mode. 
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Formally, the initial conditions for this system of equa- 
tions require some care to handle an apparent singular- 
ity, but in practice this should not be a problem. With 
a finite initial ■0(0) in Eq. @), the initial conditions for 
the Markovian closure equations are 0(0) — 0, ^(0) = 0, 
©3(0) = 0, and C(0) = Ci. For short times, we then have 
Cit) w Ci, e{t) « Cy^t. If 7JC is finite, then for short 
times we also have Mt) = Cft and 63 = {CiCf)'^/'^t'^ /2. 
It follows from Eq. (^ that 77c = - Re 77 tC//(2Ci) 
for short times, which is a consistent solution that is fi- 
nite and continuous, resolving the 0/0 ambiguity in the 
last term of Eq. (|3|). In a numerical code, it is conve- 
nient to use the initial conditions 8(0) = 0, 83(0) — 
(thus assuming the initial noise C/ — 0), C(0) = Ci, and 
A{0) — CiAt, where At is a time step smaller than any 
other relevant time scales in the problem. 



IV. FORMULATION OF THE FULL NONLINEAR 
PROBLEM AND STATISTICAL CLOSURES 

In this section we provide background on the general 
form of the nonlinear problem we are considering and 
on the general theory of statistical closures. In partic- 
ular we will write down Kraichnan's direct-interaction 
approximation, which is the starting point of our calculcu 
tion. This section borrows heavily from the BKO papertl 
(including some of their wording), but is provided for 
completeness to define our starting point. 



A. The fundamental nonlinear stochastic process 

Consider a quadratically nonlinear equation, written 
in Fourier space, for some variable "tpk'- 

d 
dt 



jf^+^k)Mt)^^ E Mi^p^rpitWgit). (35) 

k+p+q=0 



Here the time-independent coefficients of linear "damp- 
ing" Vk and mode-coupling M^pg may be complex. Given 
random initial conditions, we seek ensemble-averaged 
(or, if the system is ergodic, time-averaged) moments of 
ipkit), taking for simplicity the mean value of ipk to be 
zero. 

Many important nonlinear problems can be repre- 
sented in this form with a simple quadratic nonlinearity. 
For example, the two-dimensional Navier-Stokes equa- 
tion for neutral fluid turbulence can be written in this 
form, where ip represents the stream function such that 
the velocity v = zx V^p, and M^pq — i-pxq{q^ —p^) /k"^ . 
Other examples include Charney's barotropic vorticity 
equation for planetary fluid flow, and a class of two- 
dimensional plasma drift wave turbulence problems (such 
as the Hasegawa-Mima equation or the Terry-Horton 
equation) . Some three-dimensional one- field plasma tur- 
bulence problems can also be written in this form since 



the dominant E x B nonlinearity acts only in two dimen- 
sions perpendicular to the magnetic field. The three- 
dimensional Navier-Stokes equations and general multi- 
field plasma turbulence equations can also be written in 
the form of Eq. (js^) if ipk is considered as a vector and 
Vk and Mkpq become matrices or tensors. In fact, BKOEl 
consider covariant multiple-field formulations of the DIA 
and Markovian closures. Here we will focus on the one- 
field case, where ipk is a scalar amplitude for mode fc. 

For each fc in Eq. (|35|), the summation on the right- 
hand-side involves a sum over all possible p and q that 
satisfy the three- wave interaction fc -I- p -I- q = (this is 
sometimes expressed as fc = p-l- q, but the reality condi- 
tions tp^k = i^k has been used to rearrange it). Without 
any loss of generality one may assume the symmetry 



Mkpq — Ml 



kqp • 



(36) 



Another important symmetry possessed by many such 
systems is 



<^k.Mkpq + (JpMpqk + UqMqkp — 



(37) 



for some time- independent nonrandom real quantity Uk- 
[See Refs. ^ and 33 for the relation between this sym- 
metry and the Manley-Rowe relations for wave actions.] 
Equation (|3^) is easily shown to imply that the nonlin- 
ear terms of Eq. (|35| ) conserve the ensemble-averaged 
total generalized energy E = \'}2ik'^k{\'^Pk{t)\ )■ [The 
nonlinear terms also conserve the generalized energy in 
each individual realization, although we will be focus- 
ing on ensemble-averaged quantities, where (. . .) denotes 
ensemble-averaging.] For some problems, Eq. ( ^ ) may 
be satisfied by more than one choice of Uk] this implies 
the existence of more than one nonlinear invariant. For 
example, in the case of two-dimensional hydrodynamics, 
Eq. ( p7| ) is satisfied for both ak = and ak = k*, which 
correspond to the conservation of energy and enstrophy, 
respectively. 

We define the two-time correlation function Ck{t, t') = 
(V-'fc(^)V'fc(^')) ^iid the equal-time correlation function 
Ck{t) = Ck{t,t) (note that the two functions are distin- 
guished only by the number of arguments), so that E — 
i^^(TfcCfe(i). In stationary turbulence, the two-time 
correlation function depends on only the difference of its 
time arguments: Ck(t,t') = Ck{t — t'). The renormalized 
infinitesimal response function (nonlinear Green's func- 
tion) Rk{t,t') is the ensemble-averaged infinitesimal re- 
sponse to a source function Sk{t) added to the right-hand 
side of Eq. (^5|) for mode fc alone. As a functional deriva- 
tive. 



Rk{t,t') 



I Hk{t) 
\SSk{t') 



(38) 



We adopt the convention that the equal-time response 
function Rk{t,t) evaluates to 1/2 [although lime^o-i- 
Rk{t + e,t) = 1]. 
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B. Statistical closures; the direct-interaction 
approximation 

The starting point of our derivation will be the equa- 
tions of Kraichnan's direct-interaction approximation 
(DIA), as given in Eqs. (6-7) of BKO,El and reproduced 
below as Eqs. <^^. 

The general form of a statistical closure in the absence 
of mean fields is 

^ + i^k^ Ckit, t') + ^ Efc(t, l)Ck{i, t') 



diTk{t,t)Rl(t' ,t), (39a) 



dt 



Vk]Rk{t,t')+ / dtT.u(t,t)Rk{i,t') 



5{t~t'). 



(39b) 



While these equations (with the expressions for Sj, and 
Tk given below) are an approximate statistical solution 
to Eq. (^5|), they are the exact statistical solution to a 
generalized Langevin equation 

^ + Mt) + dt^kit, t}i;k{t) = fk{t), (40) 

where Efe is the kernel of a non-local damping/propaga- 
tion operator, and fkit,t) = {fk{t)f^{t)). These equa- 
tions specify an initial-value problem for which i is 
the initial time. 

The original nonlinearity in Eq. ( ^5|) gives rise to two 
types of terms in Eqs. (|39|): those describing nonlinear 
damping (Sfe) and one modeling nonlinear noise {^k)- 
The nonlinear damping and noise in Eqs. ( |39| ) are deter- 
mined on the basis of fully nonlinear statistics. 

The direct-interaction approximation provides specific 
approximate forms for Efc and J-'k'. 

Efc(t,i) = - J2 MkpqM;^kR;it,t)C;{t,F), (41a) 

k+p+q=0 



Tk{t.t)^\ E \Mkpq\^Cl{t,T)Cl{t,t). (41b) 

k+p+q=0 

These renormalized forms can be obtained from the for- 
mal perturbation series by retaining only selected terms. 
While there are infinitely many-ways of obtaining a renor- 
malized expression, KraichnarO has shown that most of 
the resulting closed systems of equations lead to phys- 
ically unacceptable solutions. For example, they might 
predict the physically impossible situation of a negative 
value for Ck{t,t) (i.e., a negative energy)! Such behavior 
cannot occur in the DIA or other realizable closures. 

The DIA also conserves all of the same generalized en- 
ergies (5 J2k <^k\'4'k{t)f ) that are conserved by the prim- 
itive dynamics. To show this important property, it is 



useful to write the equal-time covariance equation in the 
form 



d_ 

di 



Ck{t) + 2ReNk{t) = 2Re Fk{t), 



(42a) 



where 



Nk{t) ^ ykCk{t) - MkpqM;gk%qkit), (42b) 

k+p+q=0 



Fk{t) 



\ E 

k+p+q= 



\MkpqVQlpq{t), 



Qkpq{t) = / diRk{t, t) Cp{t, t) Cqit, t), 

to 



(42c) 



(42d) 



given initial conditions at the time t = U) (unless other^ 
wise stated, we will take = 0). As shown in B KO,E J 
the symmetries (|6|) and ( |37| ) ensure that Eq. ( f42a| ) 
conserves all quadratic nonlinear invariants of the form 
E = \^k'^kCk[t) in the dissipationless case where 
Rei/fc = 0. The Markovian closures that BKOEl devel- 
oped, and that we extend here, preserve the structure of 
Eqs. (42) and so have all of the same quadratic nonlinear 
conservation properties as the original equations. [One 
can show t hat Fk is always real, so the Re operation on 
Fk in Eq. (42a) is redundant.] 

The DIA equations (|3^) and ( plj ) provide a closed set 
of equations, but are fairly complicated because they in- 
volve convolutions over two-time functions. Their general 
numerical solution requires 0{Nf) operations, or QiN'^) 
operation&jin steady state. As described in BKCu and 
Krommes,a a Markovian approximation seeks to simplify 
this complexity by parameterizing the two-time functions 
in terms of a single decorrelation rate. Our approach here 
is essentially to generalize this to allow several rate pa- 
rameters to be used, to allow the decorrelation rate for 
Ck{t, t') to differ from the decay rate for Rk{t, t'). 



V. RESPONSE FUNCTIONS IN A STATISTICAL 
STEADY STATE 

Markovian models provide approximations that can 
simplify the integrals in Eqs. (|39|). For insight, we will 
first investigate the long-time limit where a statistical 
steady-state should be reached, so that the two-time cor- 
relation function C{t,t') and response function R{t,t') 
can depend only on the time difference t—t'. In a stalisti- 
cal steady state, all of the Markovian models in BKOl use 
a simple exponential behavior for Ck{t,t') and Rk{t,t'). 
Here we will assume the model forms 



i?mod,fc(i, t') = exp(-77fe(i - t'))H{t - t') 



(43) 



and 
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Dd,fe 



(t,t') 



Cofc exp(-77cfc(i - t')) for ^ > t\ 
Cok exp(-7]J,^(i - t')) for t < t' . 



(44) 



Note that rjk is the decay rate for the infinitesimal re- 
sponse function while r]ck is the decorrelation rate 
forCfc(t,t')- 



Inserting Eq. (41a) into Eq. ( 39b ) and using the ex- 
ponential forms of Eq. (^) and Eq. (44) in the integrals 
yields 

d 



Ol + ^k ]RkM = 5[t-t') 



J2 Mkp^M;^kCoqH{t - t') 

k+p+q=0 



X dt exp(-(7y; -I- 77jq)(t - ?) ~ Tjkit- t')). (45) 
Evaluating the integral gives 
(^^^ + '^k^Rk{t,t') = 5it-t') 

X [exp(-,/fc(t - t')) ~ exp{~{rj; + r;J,)(t - t'))] . (46) 
The solution to this equation for t > i' is 
Rk{t,t') =exp(-iykit-t')) 



E 



MkpqM*^^Coq 



k+p+q=0 'P 



Vk 



exp{-i^k{t - t')) ~ exp(-77fc(i - t')) 



Vk - Vk 

exp(-^fc(t - ^0) - exp(^(7?; + 7?g,q)(t - t')) 
Tl*p + Vhq - ^k 



(47) 



Clearly this is not strictly consistent with the simple ex- 
ponential form for Rk assun ied i n Eq. (^3|) and used to 
evaluate the integrals in Eq. ( |39b| ) . We will instead fit the 
model Eq. (^) to Eq. (p7|), in the same way that we did 
in the Langevin case for Eq. (^9|). Requiring that both 
Eq. (^) and the full Eq. (|4^) give the same weighted 
average over time (where R'^^^ f, is used as the weight to 
ensure invariance to frequency shifts) gives 



1 



Vk + Vk 



dtR:,,^,^{t,t')Rkit,t'). 



(48) 



Inserting Eq. (^) on the right-hand side, and carrying 
out a few lines of algebra, the result is 



1 



1 



Vk + Vk ^k + Vk 



MkpqMpgkCoq 



(49) 



A little rearranging leads to 



Vk = Vk 



E 

fe+p+q=0 



V*k +Vp+ Vhn ' 



(50) 



Note that this has a similar form to the steady-state de- 
cay ratftjn the DIA-based EDQNM, such as in Eq. (39b) 
of BKOtl (but with their 77* replaced by Vcq)- 

One can go through a similar calculation of Ck{t,t'), 
and calculate its weighted time average to determine the 
decorrelation rate vck- We will not do so now, as one 
can instead just take the steady-state limit of the results 
in the next section. Eq. (^) can also be obtained from 
the steady-state limit of the results in the next section, 
and so provides a useful cross-check. 

We note that there is some flexibility in the choice of 
weighting in Eq. (^). We could use C^^^ ^(t, t') as the 
weight instead of -Rmod fc(^' Either choice preserves 
Galilean invariance. Using this alternate weight, Eq. ( |4^ ) 
becomes 



kO 



Vk + Vck 



dtC*^,,^{t,t')Rk{t,t') 



(51) 



and the resulting expression for Vk is like Eq. ( pOj ) but 
with Vk on the right-hand side of Eq. (|5^) replaced by 
Vck^ which would automatically agree with the steady- 
state Vk to be defined in Eq. ([Zl]). But it turns out that 



the main steady-state results of Sec. (VII) hold with ei- 
ther choice of weights, and it seems more symmetric and 
makes more sense as a standard fitting procedure to use 
-^mod k ^ f weight for integrating Rk in Eq. ( p8[ ) . This 
raises the question of whether to use Cj^^^ ^ or R'^^^^ as 
the weight function for time averages of Ck{t,t'), as we 
will do in the next section. We can resolve this ambiguity 
by go ing b ack to the steady-state Langevin problem of 
Sec. (Ill A| ) . If one tries to use R*{t,t') as the weight in 
Eq. (|19|), so that it becomes 



Co 



Vc + V* 



dt' exp(-77*(t - t'))Cit, t'), (52) 



then one can go through the same steps used to derive 
Eq. (22) and find that in the limit of real coefficients 
it gives Vc = VVf/i'^V + Vf)- Iri the red-noise limit 
Vf <C 7?, this gives vc = '7//2, which is a factor of 2 
off from the correct result (vc = Vf) for the red noise 
limit. Thus, we will use C^^^ ^{t,t') as the weight for 
taking time-averages of Ck{t, t') and use i?*„od.fc(^' for 
time-averaging Rk{t, t'). The weighting choices might be 
reconsidered in a multi-field generalization of a Marko- 
vian closure, where the requirement of covariance may 
impose constraints on the choice of the weight functions, 
but it seems that the symmetric choices made here are 
most likely to generalize well. 
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VI. TIME-DEPENDENT MULTIPLE-RATE 
MARKOVIAN CLOSURE 

Applying these techniques in a straightforward way to 
the time-dependent DIA equations leads to the Multiple- 
Rate Markovian Closure (MRMC) equations. The two- 
time correlation function is modeled with the realizable 
form 



(53) 

(for t > t', with G„od,fc(i,t') = C*^od,ki^'^t) for t < 
and the response function is modeled as 



dirjkit}] H{t-t'). (54) 



-Rmod,fc(i, i') = exp 



Denoting Okpajt) — Qkva {t)C l^^ (1)0^^ (t), and insert- 
ing Eqs. ( |53[]54| ) into Eq. ( 42d ), we can write the equal- 
time DIA covariance equations of Eq. (p2|) as 



d_ 
dt 



Ckit) + 2Ref]kit)Ck{t) = 2Fk{t), 



(55a) 



k+p+q=0 

(55b) 

Ffc-i \Mkpq\^eipq{t)C'^/\t)Cl/\t), (55c) 

k+p+q=0 



d_ 

dt 



Qkpq + iVk + VCp + VCq)Okpq - C^/^ (t ) ^ (^ ) , 



QkpqiO) = 0. 



(55d) 



(55e) 



This is very similar to the Bowman-Krommes-Otta- 
viani Realizable Markovian Closure (RMC) (as given by 
Eqs. (66a-e) of BKOld), but with the replacement of the 
single decay/decorrelation rate of the RMC with three 
different rates in these equations. [Other Markovian 
models, such as the EDQNM closure, also use a single 
decorrelation rate parameter.] If in Eq. ( |55d| ) we replace 
Vk = fjk, rjcp = "PiVp): and rjcq = Vifjq), then these 
equations become identical to the RMC. 

To summarize the three different rates used here: 

• fjk is the nonlinear energy damping rate for the 
wave energy e qua tion for the equal-time c ovar iance 
Ckit) in Eq. (55a), and is defined in Eq. (|55b|); 



• r]k is the decay rate for the infinitesimal response 
function Rk{t,t') in Eq. (|5^), and is defined in 
Eq. (|l|); 

• and 7]ck is the decorrelation rate for Ck{t,t') in 
Eq. (p3|), and is defined in Eq. (|69|). 

To determine rjk{t) and rjck{t), we follow a similar 
procedure as we did fo r the time-dependent Langevin 
equation in Sec. ( Ill C| ). Define Ak{t) as the following 
weighted time-average of Rk 

Ak{t)^ I dt'Rl^^^^^{t,t')Rk{t,t'). (56) 

If Rk{t, t') = i?mod,fc(i, t') as given by Eq. (H), then 
dAk 



dt 



1 - {'nk+Vk)Ak 



(57) 



while ii Rk{t,t') satisfies Eqs.( p9b|j4Ta| ), then 
dAk 



dt 
where 



^-Wk + Vk)Ak 

+ ^ ^^kpqMpq^Cq^^Ql pqk, 
k+p+q=0 



C*Jt,t) 



(58) 



' dt' R*^,^^kit,t') / dt^i}^R;{t,t)Rk{t,t'). (59) 
Jt' C'l'^it) 



It is often more convenient to work with the differential 



version of this, which, after using Eqs. (53-5J) to replace 
Cq{t,t') and Rp{t,t') with their model forms, is 



de 



l.pqk 



dt 



= -ivl + vhq + v;)Qlpqk + c'q^'Mt) (60) 



(with the initial condition 6i.fcpq(0) = 0)). Requiring 
that Eq. (57) and Eq. (58) be equivalent determines ijk 
to be 



Vk 



= - 7- E MkpqM;^^cl/^ei 



pqk ■ 



(61) 



k+p+q=0 



The calculation of rjck proceeds in a similar way. 
Ack{t) is defined as a weighted time integral of Ck{t, t'): 



Ack{t)= / dt' C:^,,^{t,t')Ck{t.t'). 



(62) 



If Ck{t,t') in this integral is replaced by Cinod,k{t,t') as 
given by Eq. (p3|), then 



dt 



1 dCkjt) 
Ck{t) dt 



A, 



Ck 



iVck + Vck)A. 



Ck, 



(63) 
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(where we make the time dependence of Ck (t) explicit to 
distinguish it from the two-time Ckjt, t')) . If the exact 
dynamics for Ck{t,t') given by Eqs. (39a, 4l|) are used, 
then 



Ck 



dt 



cl{t) 



dCkit) 



2Ck(t) dt 



Ack - ivhk + '^k)A, 



E 

k+p+q=0 



A,f, i\T* r'^/^r'i/^c)* 



\Mkpq\^Cl''cy^Cl'-'Qlkpq, (64) 



fc+p+q=0 



where 



e 



At) 



10 

and 



^t :;tjt. 



dt 



C*{t,t) 

1/2 



Ck' (t) Jo Cq'^it) 



R*Jt,t)Ck{lt') (65) 



e 



At) 



3,kpq ' 

c'J'it) 







^.c*it,t)c;{t,i) 



RUt\t). (66) 



ci^^{t)c^/\ty 

Using Eqs. ( p3[]5^ ), the differential versions of these are 



2,pgfc 



dt 



-ivhk + vhq + v;)Q 



2,pqk 



'Ck{t)e;,kit) 



Cq^'jt) 

cl^\t) 



Ackit) (67) 



and 



3,kpq 



dt 



^iVCk + VCq + VCp)'^3.,kpq 

+ci^'mipqit). 



(68) 



The quantity rjck is then determined by equating 
Eq. (|63|) and Eq. (p4|), yielding 



VCk = Vk 

1 



1 dCk{t) 



Ack 



2Ck{t) dt 

^kpq^^J^pqk'-^k 



J2 Mk„.ML^c'J^cl^^e* 



2,pqk 



fc+p+q=0 



2Ack 



E IMkpql'cl^'c^'cy'ei 



3,fcpq ' 



(69) 



k+p+q=0 



where Eq. (55a) could be used to eliminate dCk{t)/dt. As 
in Eq. ( p3| ) for the case of the time-dependent Langevin 
equation, while there are effects in this equation that will 
tend to give Tlerjck > 0, it may be necessary to modify 
this equation to enforce realizability in all cases. This is 
done by replacing this equation, of the form rjck — RHS, 
with r]ck = 7'(RHS). Note that it is only Re r]ck > that 



is needed for realizability, while Re?7fc can transiently 
go negative (as it does in two-dimensional hydrodynam- 
ics because of the inverse cascade, or in some plasma 
problem s whpr e the zonal flows may become nonlineasly 
unstablcEiyO). This is simflar to BKO's treatment Jd 

The complete set of equations that constitutes the 
Multiple-Rate Markovian Closure (MRMC) are Eqs. (|5|) 
for the equal-time covariance Ck{t) and related quan- 
tities, Eqs. ( [57|j60| , |6l| ) for quantities related to the re- 
sponse function, and Eqs. ( |63| , |67[]6^ ) for quantities re- 
lated to the two-time correlation function. The MRMC 
extends the RMC to make less restrictive assumptions 
and include additional effects, but at the expense of a 
few new parameters. In addition to replacing the sin- 
gle decay/decorrelation rate of the RMC with 3 different 
rates, rjk, fjk, and rjck, it also replaces the single triad 
interaction time of the RMC with 4 different triad inter- 
action times, Qkpq, Qi,kpq, 62,fcpq, and Qs^kpq- Each 
of these triad interaction times has a different weighting 
of response functions and two-time correlation functions. 
While this increases the complexity some, the overall 
computational scaling of this system is still 0{Nt), a 
significant improvement over the 0{Nf) or 0{Nf) scal- 
ing of the DIA. 



VII. PROPERTIES OF THE MULTIPLE-RATE 
MARKOVIAN CLOSURE 

In a steady-state limit, Eq. (|6l|) simplifies to 



Vk = Vk 



E 



fe+p+<j=0 



MkpqM;^^c^ 

V*k+Vp + Vhq ' 



(70) 



The steady-state balance RerjkCk = Fk from Eq. (55a) 
simplifies to 



Ck Rc Vk 



k+p+q=0 ^Ck +Vp + Vhq 
\Mkpq\^CrjC, 



E * * * 

fc+p+q=0 ^fc + ^Cp + VCq 



(71) 



[Note the subtle notational differences: the expression 
for r]k becomes the expression for fjk if rj^ on the RHS of 
Eq. ( [70| ) is replaced by 7]^^ .] Finally, Eq. ( |69| ) reduces to 



Vck = Vk- ivck + vhk) 

E 



1 

2Cfc 



E 

k+p+q=Q 
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Mkpq^^pqkCq 



Ck 



Vcn 



kpq\ CpCq 



^0 (V*ck+V*cp + V*cq)il*k+V*cp + V*cq) 



(72) 



Thus the decorrelation rate -qck equals the response func- 
tion decay rate rjk plus the two correction terms in brack- 
ets. For the simple steady-state non-wave case with real 
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and positive ry's, the second correction term will cause 
Tjck to decrease (as expected for non- white noise), while 
the first term will usually have an offsetting opposite sign 
and cause rjck to increase. The origin of th ese two terms 
can be traced back to the DIA Eq. (39a). The second 
correction term corresponds to the usual effects of non- 
white noise (related to the integral involving J-k (t, t) in 
Eq. (39a)), but the first correction term in Eq. ( |7^ ) is 
related to the time-history integra l inv olving the renor- 
malized propagator Efc(t, t) in Eq. (39a). Thus non- white 
fluctuations in other modes Cq{t, i) not only change the 
noise term for the k mode, but also change the effec- 
tive damping from the time-history integral, broaden- 
ing the width of Sfc(f, i) in time (if the fluc tuations Cg 
were treated as white noise, then Eq. (41a) would give 
Sfc(i,F) cx S{t-t)). 

An important property to demonstrate is that in ther- 
mal equilibrium it is possible for these two terms to can- 
cel exactly. Then the decorrelation rate and response 
function decay rate are equivalent, r]ck — Vk, and the 
fluctuation-dissipation theorem is satisfied. To demon- 
strate that this is true, we assume the result (r?c/c ~ Vk) 
to simplify some of the equations and then show that this 
is a self-consistent assumption. (Note also that if rjck = 
i]k, then fjk — rjk also.) Splitting the first summation in 
brackets in Eq. (|7^) into two equal parts and interchang- 
ing the p and q labels for one of these parts (i.e., using 
an identity of the form SGfcpg = YuGkpql^ + '^Gkqpl^), 
the terms in brackets in Eq. (|7^) can be written as 



— y — 



M, 



kpq 



X {Mpg^CqCk + M*p,^CpCk + MlpqCpCq). (73) 

In thermal equilibrium, the spectrum Ck is given by 
equipartition among modes of a generalized energy-like 
conserved quantity. Consider an equipartition spectrum 
of the form Ck — 1/Afc, where Afc = a(*V^'\ the ct^*'' 
are the coefficients in Eq. ( |37| ) (related to the quadratic 
invariants), and a'^^^ are determined by the initial con- 
ditions. Substituting Ck = 1/Afe, Cp = l/Ap , and 



Cq = 1/Aq into E q. (|73| ), and using Eqs. (|36[-|37|), one 



can show that Eq. (|73|) indeed vanishes, so that Eq. (|72 
simplifies to rjck — Vk- The proof that Ck = 1/Afc is a 
solution of the steady-state Eq. (|7l| ) proceeds in a simi- 
lar way, interchanging the p and q labels for half of the 
summation on the left-hand side of Eq. (^l]), and not- 
ing that Rei/fc = in an isolated thermal system, etc. 
Rigorously, this only shows that the equipartition spec- 
trum Ck = 1/Afc is an equilibrium solution. This paper 
doesn't demonstrate that it is a stable equilibrium or 
that mixing dynamics will necessarily relax to this state. 
For a discussion of the Gibbs-type H theorem that leads 
to this result, see Appendix H of Ref. ^ and Refs. ^ 
and 1^. It is significant to note that the thermal equilib- 
rium result holds even if the number of modes is small, 
and it does not assume that the noise spectrum is white. 



This is unlike a simple Langevin equation of the form 
of Eq. ^ (which has a local damping term in contrast 
to the time- history integral of Eq. (^0|)), where the two- 
time correlation function and the infinitesimal response 
function are proportional only if the noise is white. 

We next estimate the importance of the correction 
terms in the decorrelation rate for an inertial range of 
a turbulent steady state, such as in two-dimensional hy- 
drodynamics. Typically most of the energy is at long 
wavelengths {Cq is peaked at sufficiently low q), so that 
the dominant contributions to the sums in Eq. ( [7^ ) come 
from long wavelengths: \q\ <^ \k\ in the first sum, 
and \q\ ^ |fc| or |p| — |fc -I- q| ^ |fc| in the second 
sum. This means that one can approximate the de- 
nominators in the sums of Eq. (|7^) using, for example, 
ivck + Vp + VCq) « ivck+Vk) (since r]k and r]ck are typ- 
ically increasing functions of k). Similar approximations 
give {fjk - Vk) ~ {Vk- '^k)2ril/^{ri*(j^ + vD- Usmg the 
steady-state relation Fk — RefjkCk from Eq. ( [7l| ) and 
the disparate scale approximations to rewrite the second 
sum in Eq. (^2|) in terms of fjk , and allowing finite dissi- 
pation but ignoring wave dynamics (so that Vk and the 
various 77 coefficients are real) , one can show that Eq. ( [T^ ) 
simplifies in this disparate scale limit to 



VCk = Vk - Vk 



'2r]k{rik - Vk){rik - Vck) 
{rick+VkY 



(74) 



This gives a cubic equation for r/cfc ■ For i^k — the roots 
are T]ck = Vk and vck = (~1 ± \/2)?7fc- Our speculation 
is that rjck — Vk will be the usual case in a steady- 
state inertial range. (This appears reasonable, but it 
might require numerical simulations to test it more defini- 
tively.) The other roots are probably unstable equilibria, 
so that any perturbation away from it would eventually 
approach the stable root, or may only be relevant in tran- 
sient inverse-cascade cases where Re 77fc < (Re vck > 
being required to satisfy realizability) . 

Thus the two correction terms in Eq. (|7^) again exactly 
cancel each other (assuming the root choice made above), 
leading to vck — Vk and the result that non-white-noise 
corrections are asymptotically unimportant in a wide in- 
ertial range (fe large compared to the long-wavelength 
energy-containing wave number scale fco). However, this 
may be an artifact of the problem that the underlying 
DIA, on which the MRMC is based, does not s*tisf 
dom Galilean invariance. As is well knownE2IE3 
the reason the DIA predicts a slightly different spec- 
trum {E{k) - fc-3/2 in the energy cascade inertial range) 
than the Kolmogorov result [E[k) ~ k~^/'^) is because 
of this lost random Galilean invariance. [The stan- 
dard definitions for two-dimensional hydrodynamics use 
E{k) ~ A:'^C|fc| when ijjk represents the stream function, 
so that the total energy is / dkE{k), a one-dimensional 
integral over the magnitude of fe.] The magnitude of 
this discrepancy between the DIA and dimensionally self- 
similar predictions is calculated for a general equation of 
the form Eq. (^5|) in Appendix |^. 
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The underlying reason for this failure of the DIA is 
that the nonlinear damping and noise terms (the left- and 
right-hand sides of Eq. (71)) are dominated by contribu- 
tions from the energy at long wavelengths. A random- 
Galilean invariant theory should depend only on the 
shear of longer wavelength modes (as rjk does in Orszag's 
phenomenological EQDNM) and the most energetically 
significant interactions should occur among comparable 
scales (|q| ~ \p\ ~ |fc|). Then the disparate scale ap- 
proximations that led to Eq. ( |7^ ) would no longer be 
valid. In such a case, it would seem unlikely that the two 
terms in Eq. (72) would still exactly cancel, and there 
would probably be some difference between the decorre- 
lation rate rjck and the decay rate rjk- It would therefore 
be interesting to try to apply the techniques developed 
here (for allowing multiple rates) to other starting equa- 
tions that respect random Galilean invariance, such as 
the Lagrangian-history DIA, test-field model, or renor- 
malization group methods. 

A regime where the correction terms might not can- 
cel each other, and the differences between rjck and rjk 
might be significant, even with the DIA's overemphasis 
of long-wavelength contributions to the eddy turnover 
rate, is in ITG/drift-wave plasma turbulence, where the 
spectrum can often be anisotropic and have strong wave 
effects. That is, Vk can be complex, with unstable modes 
in some directions and damped modes in others, so that 
ffc and C'k vary strongly with the direction of k. Some 
plasma cases have a reduced range of relevant nonlin- 
early interacting scales, and the simplifications of dis- 
parate scales in an inertial range used to derive Eq. ( |7^ ) 
are not appropriate. The corrections might also be im- 
portant in non-steady-state transient cases (such as zonal 
flows with predator-prey dynamics) or in other regimes 
where interactions between comparable scales dominate. 
Evaluating the difference between the decorrelation rate 
Vck and the decay rate r]k in more general cases such as 
these probably requires a numerical treatment. 

Finally, it is useful to demonstrate that the Multiple- 
Rate Markovian Closure approximation preserves realiz- 
ability, which turns out to require one additional con- 
straint. The MRMC equations (|5^) have the underlying 
Langevin equation 



we find that if the two-time statistics of fk satisfy 



dt 



(75) 



where fjk is given by Eq. ( |55b| ). The statistics 
that fk must satisfy can be found by comparing 
the solution for such a Langevin equation, given 
by Eq. (||), with Eqs. (|55|), finding the constraint 
ReJ*dtR*{t,t)C}{t,t) = Fk, where Fk is given by 



Eq. (55c) and R{t,t) = cxp(— dt" f]k{t")) is the prop- 
agator for Eq. ([75|) . Using an integral form for Qkpq 
(similar to Eq. ([42d|)), 



Qkpq{t) 



iT r, /j- T\ ^mod.pit,t) Caiod.qit^t) 

at Kmod,k[t, t) ■ 



Cf{t, t) — exp 



dt [f-jkit ) - 7]k{t )) 



Cl'\t)Cq''{t) 



1/2, 



X— ^ ^ \Mkpq'\^CTjiod,p{t,t)Cijiod,qit,t) 
k+p+q=0 

(for t > f), then the MRMC is the statistical solution 
of Eq. ([j'5|). As shown in Theorem 1 of Appendix B 
of BKOEr(and as can be inferred from considering the 
statistics of f{t) = g{t)h{t), where g and h are sta- 
tistically independent), a product of realizable correla- 
tion functions is also a realizable correlation function. 
Cmod.pitjt') and Cniod,q(i, i') are individually realizable 
because Rerjck > for all k. So in order to guarantee 
realizability of Cf{t,t), we need to impose the additional 
condition that He rjk > Tie fjk- This constraint seems 
physically reasonable. The parameter r/j, measures the 
decay rate for the ensemble averaged response {5ipk{t)), 
which can decay either as energy is nonlinear transferred 
out of mode k or as the energy that is in Sij^k becomes 
randomly phased. The quantity fjk used in Eq. ( ^5| ) mea- 
sures only the rate at which net energy (regardless of 
phase) is transferred out of mode k into other modes, 
so it would seem reasonable that fjk < ?7fc will naturally 
result. 



VIII. CONCLUSIONS 

In summary, we have demonstrated a method for ex- 
tending Markovian approximations of the DIA, to allow 
the decorrelation rate for fluctuations to differ from the 
decay rate for the infinitesimal response function (the 
renormalized Green's function or nonlinear propagator). 
This can give a more accurate treatment of various ef- 
fects such as non-white-noise forcing terms. In practice, 
the corrections to the decorrelation rate are modest, at 
least in isotropic non-wave cases, since the decorrelation 
rate of the noise is usually comparable to, if not much 
larger than, the decay rate for the response function. 
For example, if rjf = 77 in the simple Langevin example 
of Eq. (p2[), then the decorrelation rate is ~ 60% lower 
than its white-noise value. Furthermore, the Multiple- 
Rate Markovian Closure Eq. (|7^) for the full DIA con- 
tains an offsetting term that can increase rjck, so the net 
result is less clear. This is because the DIA is related 
to a generalized Langevin equation Eq. (|40|), where non- 
white fluctuations modify not only the noise term (which 
tends to reduce the decorrelation rate) but also modify 
the renormalized propagator in the time-history integral 
(which tends to increase the decorrelation rate). 

We have demonstrated that these two terms in fact ex- 
actly cancel each other as they should in thermal equilib- 
rium where the fluctuation-dissipation theorem applies. 
We have also found another case, that of a wide inertial 
range with no waves, where it is possible for these two 
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corrections to offset each other exactly, so that the decor- 
relation rate and the decay rates become equal. How- 
ever, this may be an artifact of the loss of Galilean in- 
variance in the Eulerian DIA, where modes in the iner- 
tial range nonlincarly interact predominantly with long 
wavelength modes. Thus it would be interesting to try 
to apply the techniques developed in this paper to other 
renormalized statistical theories, in which the dominant 
nonlinear interactions in an inertial range are between 
comparable scales instead of disparate scales and which 
properly reproduce Kolmogorov's E{k) oc k~^^^ inertial- 
range energy spectrum instead of the Eulerian DIA's 
E{k) oc Single-rate Markovian approximations 

have been applied in the past to other renormalized sta- 
tistical theoriealj and white-noise assumptions have also 
been emplcaied in renormalization group calculations of 
turbulence.tj An interesting question is whether there is 
some way to generalize such calculations to allow for mul- 
tiple rates and non- white noise as considered here. An- 
other question is whether multiple-rate extensions might 
modify subgrid turbulence models. [Such corrections 
would probably be important only at short scales near 
the transition from resolved to unresolved scales.] 

Even in the context of an Eulerian DIA-based theory, 
there may be some regimes where the multiple-rate cor- 
rections in this paper may be important and warrant 
further investigation. These might include cases where 
non-steady-state dynamics are important (i.e., predator- 
prey oscillations between different parts of the spectrum, 
such as between drift waves and zonal flows), or where 
interactions between comparable |fc| scales are more im- 
portant, such as might occur in anisotropic plasma tur- 
bulence with wave dynamics and with instability growth 
rates or Landau damping rates that vary strongly with 
the magnitude and direction of the wavenumber. One 
could test whether these corrections are important or 
negligible iaryarious regimes by looking at 3-mode cou- 
pling caseSjtrB or by numerically comparing with the DIA 
or direct numerical simulations. 

The complete set of equations that constitutes the 
Multiple-Rate Markovian Closure (MRMC) are summa- 
rized in the final paragraph of Sec. (VI). The MRMC ex^ 
tends the Realizable Markovian Closure (RMC) of BKOEI 
to allow various nonlinear rates and interaction times to 
differ. The single decay /decorrelation rate of the RMC is 
replaced with 3 different rates, ijk (the response function 
decay rate) , r]ck (the decorrelation rate for the two-time 
correlation function), and fjk (the energy damping rate). 
The triad interaction time of the RMC is replaced with 4 
different triad interaction times with various weightings 
of decorrelation and decay rates. While this increases the 
complexity of the equations somewhat, the main compu- 
tational advantages of a local-in-time Markovian closure 
relative to the non-local-in-time DIA are retained. 
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APPENDIX A: REALIZABILITY OF A 
PARTICULAR TWO-POINT CORRELATION 
FUNCTION 

In theorem 2. of their Appendix B, Bowman, Krommes, 
and OttavianiEl show one way to prove that a two-point 
correlation function of the form of Eq. ( p7| ) is "realiz- 
able" (if Keric{t) > is satisfied almost everywhere). 
Realizability means that this two-point correlation func- 
tion is the exact solution to some underlying stochastic 
problem, such as a Langevin equation. In the absence of 
realizability, non-physical difficulties can sometimes de- 
velop, such as the predicted energy C{t) = C{t,t) going 
negative or diverging. Here we present an alternate proof 
that Eq. ( p?!) is realizable. 

Consider the standard Langevin equation with time- 
dependent coefRcients, but in the white-noise limit 
(/(0/*(^')> = 2D{t)6{t ~ t'). Then Eq. (|) simplifies 
to 



dCjt) 
dt 



2Re r){t)C{t) = 2D{t), 



(Al) 



while the equation for the two-time correlation function, 
Eq. (I), becomes just dC{t,t') / dt + r]{t)C{t,t') = for 
t > t', with the boundary condition C{t',t') — C{t'). 
Taking the time derivative of Eq. (|27| ) gives 



d_ 

dt 



+ f?c(t) anod(i,i') 



1 dC{t) 
2C{t) dt 



Cmod(^: t )• 

(A2) 



If C{t,t') = Ctjiod{t,t'), then these last two equations 
give V ^ Vc - {dC{t)/dt)/{2C{t)). Using Eq. (^), 
this becomes r]c{t) = r]{t) - Rer]{t) + D{t)/C{t). It 
is interesting to note that this ensures Re 77c ^ even 
if Re 77 < 0. These equations can be rearranged to 
give lmrj{t) = Imrjcit), D{t) = C{t)Rericit), and 
Rer]{t) = Rcr^c{t) - {d\ogC{t) / dt) /2. Thus, given any 
3 arbitrary functions C{t) > 0, Re 77(7 (t), and lmrjc{t) 
that determine the model Eq. (p7|), it is possible to find 
a white-noise Langevin equation for which it is the exact 
solution (as long as Re rye > so that D > Q). Con- 
versely, for any arbitrary complex 7]{t) and real D{t) >Q 
that specify a white-noise Langevin problem, one can find 
a corresponding solution of the form Eq. (p7|) . 



18 



It is interesting to note that C{t,t') = 
C(i') exp[— /j, dt" r]{t")] (for t > t') is also an exact so- 
lution for this same white-noise Langevin problem. This 
form is valid for arbitrary rjit) {gysu Re 77 < 0). How- 
ever, BKOH and references thereirJid indicate that this 
fails to preserve realizability when used in the context of 
Markovian approximations to the DIA, so they instead 
useEq. (H). 

On a related topic, BKOl showed that their re- 
alizable Markovian closure (RMC), as given by their 
Eqs. (66a-e), has an underlying Langevin representation 
given by their Eq. (67) with a two-time noise correla- 
tion function {fk{t)fl{t')) of the form of their Eq. (64), 
which is not necessarily white noise. However, other 
two-time noise correlation functions can also give the 
same equal-time statistics equivalent to their Eq. (66a). 
This requires Fk{t) = Re J* dt {fk{t)f^{t))Rl{t,t), where 
Fk{t) is the noise term in their Eq. (66a). For a 
case where Fk{t) is always positive, then the RMC is 
also equivalent to a Langevin representation with white- 
noise, (/fc(0/fc(i')> = ^Df (1)6(1 - t'), where Df(t) 



white and non-white noise can give the same equal-time 
equations for C(t), they will give different results for the 
two-time correlation function C(t,t'). However, there 
can be cases where Fk{t) < 0, for which a realizable 
Langevin representation must use non-white noise, as in 
their Eq. (64). [Note that while C{t) > is a funda- 
mental requirement preserved by a realizable theory, the 
"triad interaction time" Re &kpq niay go negative. An 
example, similar to Eq. (47) of BKO,El can be constructed 
for the realizable Qkpq of Eq. (66d) of BKO in the limit 
of constant Cp and Cq with rik = rjp + r/g = p + ia.] 



APPENDIX B: FITTING MODELS TO THE 
TWO-TIME CORRELATION FUNCTION 

Conceptually the process of fitting an exponential 
model of decorrelation to the actual two-time correla- 
tion function seems straightforward. But as described 
in Sec. ( Ill A ) and Sec. (0), there are various choices 



one could make in the weights used to fit the models. 
Galilean invariancc imposes some constraints, but does 
not completely constrain the problem. In this appendix 
we further describe some options and our choices. 

Consider the following measure of the error between 
the actual two-time correlation function and a model cor- 
relation function: 



S(t) = / dt'\C(t,t') 
Jo 



Cmod(^, t')\'' 



(Bl) 



We win assume Cmod(t,'t) is of the form of Eq. (pT]). The 
equal time correlation function C{t) — C{t,t) is already 
specified, so our task is to choose ric(t) in Eq. ( p7| ) in such 
a way as to minimize the squared error S. We want to 
stay in a Markovian framework, where ric(t) depends on 



parameters only from the present time. We assume that 
ric(t') for times t' <t has already been chosen optimally. 
But we can choose rjc(t) at the present time so that the 
extrapolation of S(t) into the future is minimized. That 
is, we want to minimize dS/ dt, which, after using Eq. (|^) 
for dC(t,t')/dt and Eq. ^) to evaluate dCn,odit,t')/dt, 



9S „ /■* , 



f](t)C(t,t')+ / dtR*(t',t)C}(t,t) 



1 dC(t) 



c„,od(t, t') + 'nc(t)a^od{t,t') 



2C(t) dt 
x(C*(t,0-C*.od(i:0)+c.c., 



(B2) 



where c.c. indicates the complex conjugate of the pre- 
vious expression. Separately minimizing dS/dt with re- 
spect to the real part rjcr and imaginary part rjQ^ of rjc (t) 
(i.e., set d(dS/dt)/ drier = 0, and ilvci^ d(dS / dt) / ^r]c^ = 
0) leads to the requirement that 



l^eJ2^+p+q=o\Mkpq\^&kpqCp^^Cl^^- While both 



* dt' t')a^od(t, t') = /* dt' C*^,S. t')C(t, t'). 

Jo 



(B3) 



[Note that when evaluating derivatives of Eq. (|B^) with 
respect to r]cr and ri a, i t is only the explicit appear- 
ance of r]c(t) in Eq. (|B2| ) that is important. The pa- 
rameter ric(t) also appears implicitly via the definition 
of Cmod(t, t'), but there it has an impact on the integral 
defining dS/dt only through a set of measure zero, and 
so can be neglected as long as ric(t)js bounded.] 

In the steady-state limit, Eq. (p3) is equivalen t to 
Eq. (|l^). For a time-dependent case, consider Eq. (B3) 
as providing a constraint of the form Amod(^) = ^(^)- 
Assuming that this has already been satisfied for earlier 
times, we want it to remain satisfied for future times, 
i.e., we need to require that dA^md/dt — dA/dt. This 
is precisely what we are doing when we set Eq. ( ^9|) and 
Eq. ( |30| ) to be equal, and it leads to a formula for r]c(t) 
at the present time that minimizes the errors as time 
advances. 

The same procedures as described here are used in fit- 
ting a model response function i?mod(^, t') of the form of 
Eq. (54) to the actual response function, leading to the 
constraint 

t pt 
dt' R^oditTi')Rniod(t,t') — / dt' R^^^{t,t')R{t,t') 

Jo 

(B4) 

As mentioned at the end of Sec. (^, R'^^^ in this ex- 
pression could be replaced with C^^^ and one would still 
get an expression defining 77 that was Galilean invariant. 
However, Eq. (B4) seems to make more sense as a least- 
squares best fit of i?mod to R, and that is the choice we 
have made. 
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But consider Eq. (Bl) in the steady-state limit where 
rjc is a constant and C(t, t') depends only on t — i', 



So 



dt'|C(<,i')-C'oe"''^'*"*')p 



(B5) 



It is straightforward to show that choosing r/c to mini- 
mize the total squared error Sq leads to the condition 

r dt' e-''c(*^*')Coe^'"=(*^*')(t - t') 

dt' e-''cit-i')c{t,t')[t-t'). (B6) 



Note that this differs from Eq. ( p^ ) by an additional fac- 
tor of {t — t'), which weights errors at larger time separa- 
tion more strongly. Including an extra weighting factor 
of {t — t') in Eq. ( |l9| ) might help to refine the model, 
particularly for cases such as in Fig. (5), where the short 
time behavior is reasonable but the long-time fit needs 
improvement. 

It is perhaps not surprising that optimizing a constant 
•qc to minimize the global error gives a somewhat dif- 
ferent result than optimizing rjc{t) to minimize the local 
error dS/dt. In order for the time-dependent fitting pro- 
cedures to reproduce this steady-state result, one could 
modify Eq. (Bl) by multiplying the integrand by a fac- 
tor oi {t — t'). Working through the derivation, one finds 
that the integrands in Eq. (B3) would be modified to also 
have an additional weighting factor of {t — t'). Thus one 
might be able to improve the results in this paper some 
by including an extra weighting of {t — t') in the appro- 
priate places, Eq. (||), Eq. Eq. (||), and Eq. (|62[), 
and working through the derivations to see the modified 
results. While such modifications could lead to an im- 
proved model, and would be interesting for future work, 
one should realize that the dynamics are complicated and 
no choice of weights is perfect. For example, what one 
really wants is a best fit model for the triad interaction 
times which are weighted by interactions between three 
modes as given in Eq. (42d), not necessarily best fits for 
the decorrelation rates of just individual modes. Proba- 
bly a higher priority for future work is to use a starting 
set of equations that satisfy random Galilean invariance, 
so that interactions with large scales are not overempha- 
sized as they are in the Eulerian DIA. 



APPENDIX C: INERTIAL-RANGE SCALING OF 
DIA-BASED CLOSURES 

Here we determine steady-state self-similar inertial- 
range solutions in d dimensions to closures of 
the form ( ^2[ ) in an unbounded domain (so that 
Efc+p+q=o ^ lA^dpdq = J dpdqS{k+p + q)), taking 
the initial time to — — oo. This extends previous deriva- 
tions in the literature to self-similar spectra consistent 



with generic DIA-based closures ( p^ ) of the quadrati- 
cally nonlinear equation (|35|), arising from the cascade 

of a generalized energy i crfc|V'fc(OI ■ Assuming self- 
similar scalings of the mode-coupling and statistical vari- 
ables, our derivation requires only the additional condi- 
tion ( p^ ) , which is somewhat weaker than statistical sta- 
tionarity. 

The turbulence could be forced with a linear instabil- 
ity, incorporated with dissipation into the linear coeffi- 
cient , or else a random force could be added to the 
right-hand side of Eq. ( 42a ). By definition, both the ex- 
ternal forcing and dissipation Vk vanish in the inertial 
range. The symmetry ( |37| ) then implies that the nonlin- 
ear terms in Eq. (42a), weighted by ak, must balance. It 
is convenient to define 

Sk ^ <JkRe{Fk - Nk) 

^■Re / dpdqukMkpqMkpqQkpq 



/A 

Re 



= -iRe 



dp dq <ykMkpqM*g^e*pg^ 



dpdqMkpq{(TpM*g^ + <JqM*^p)ekpq 



Re / dpdqcTkMkpqMpgk'dpqk 

'Afc 



= - Re 



dp dq apMkpqMpqk^kpq 



+ Re / dp dq cTkMkpqM;^^%^^ 

JAk 

= Re / dpdqMkpqMpqk{crkQpqk~ f^pQkpq)- 

We seek self-similar solutions of the DIA that obey the 
scalings (for A > 0) 



M 



XkjXp.Xq 



kpqj 



Rxk{t.t') ^ Rk{t,t- X-'{t-t')), 



Cxk{t,t')^X''Ck{t,t-X-\t-t')), 



(Cla) 



(Clb) 



(Clc) 



(Cld) 



so that, upon making the change of variables s ^ t — 
\-\t-^ in Eq. (|2d|). 



e 



Xk,Xp,Xq 



kpq- 



(Cle) 



Once we have determined suitable values of the scaling 
exponent n, we may compute the wavenumber exponent 
fi for the energy spectrum E{k) ~ e°'k^ . If the total 
energy E is related to the correlation function Ck of the 
fundamental variable ip hy E — J dkk'^Ck — J dkE{k), 
then P = d— 1 + j + n. 
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Following Ref |Tl|, we will use the change of variables 
z = fc^/p, w = fcg/p to determine values of the exponents 
I and n for which the angular average S(}i) of Su vanishes. 
In terms of the scaling factor A = kjz we note that k — 
Az, "p — \k, and q — Xw. Letting z = zp and w ~ wq, we 
may then express dpdq — X'^'^dz dw and S{k + p + q) = 
X^'^S{zk + kp + w). Hence, upon interchanging p and k 
in the integration, we deduce 



S{k) = 



dk Sk ^ J dk J dz dw A 



3d~d+2m+s+e+2n 



Rej dkj dzdwX^^+^"'+'+'+^" 



k.. 



= -S{k), 
provided that 

2d + 2m + s + i + 2n = 0. 



(C2) 



The condition (C2) gu aran tees that the angle-averaged 
nonlinear terms in Eq. (42a ) will balance in a steady state 
and lead to an inertial range. 

The exponent £ can be determined by integrating the 
DIA response function equation 



^Rk {t, t')- f dt f dp dq MkpqM;^^ 
xR;{t, t) C* (t, t) Rk (t, t')^6{t~t'), 
over all t', using the steady-state condition 



lim — / dt' R{t,t') 



0. 



(C3) 



(C4) 



One obtains 



dt / dpdqMkpqM*^^ 



/OO 
dt'Rk{t,t') = 1. 
-OO 



(C5) 



Upon replacing k by Xk (for any constant A) and exploit- 
ing the self-similar scalings given in Eqs. (CI), we make 
the change of variable s' = t — X~^(t — t') to obtain 



<.d+2m+t+n 



dt I dpdqMkpqMpgf, 



xi?;((X),s)C*(oo,s) 



ds'Rkit,s') = 1, 



where s = t — X^^{t — t). The integral over t is domi- 
nated by contributions from large t, for which the inte- 
gral over s' asymptotically ap proa ches a constant (with 
respect to i), according to Eq. (C4). Hence, after making 



a final change of variable s fro m t to s, we see that the 
balance expressed in Eq. (C5) is recovered if 



A 



d+2m+2i+n _ 



= 1, 



(C6) 

-{d-\-n)/2 — m. If one 
inserts this result into Eq. (C2), one obtains the Kol- 
mogorov scalings 



from which we conclude that 



1 



n = —d — -(m + s). 



/3 = 7-l- -(m + s). 



(C7a) 



(C7b) 



(C7c) 



Alternatively, one could adopt instead of Eq. (|Cj) the 
stronger condition of statistical stationarity, Rk {t, t') = 
TZk{t - t') and Ck{t,t') = Ck{t - t'). Equatio n (|cq ) is 
then readily seen to follow directly from Eq. (|C3|). In 



either case we have only shown that Eq. (C6) is a neces- 
sary condition for self-similar solutions of the form (CI) 
to exist. In order that these solutions actually satisfy 
Eq. (p3|), it is also necessary at the very least that the 
wavenumber integral in Eq. (C2) converges. 

Unfortunately, the scaling expressed in Eq. (|C6|) o f- 
ten leads to a divergence of the q integral in Eq. (|C3|), 
preventing self-similar solutions from existing. Typically, 
the mode-coupling coefhcients M^.-k-q.q asymptotically 
approach a constant as q goes to zero while k is held fixed. 
Upon performing the p integration in Eq. ( |C3| ) , we then 
see that the q integrand will scale like q'^^^q {t, t) for 
small q. If Cq asymptotically scales as g", t hen the inte- 



grand will scale like qd-^+n^ g^^^ ( C7b ) implies that 
d— l + n = — 1 — 2(m + s)/3. Normally m + s > (see 
Table ^; in these cases there would be a divergence of 
the q integrajr^]i| Eq. ( |C3| ) if self-similar solutions really 
were to exist J23'EI 

This divergence indicates that the dominant contribu- 
tions to the eddy-turnover time come from the energy 
spectrum at large scales, where self-similarity no longer 
holds. (For this reason, the DIA is not invariant to ran- 
dom Galilean transformations.) The actual value of the 
scaling £ that appears in the DIA response must be cal- 
culated by taking into account that Cq does not actually 
behave as g" for small q. The DIA equations apply to the 
case of zero mean flow, where the energy spectrum goes to 
zero at smal l wa venumbers. This means that the integra- 
tion in Eq. (C3) must be effectively cut off at some fixed 
large scale wavenumber ko. The introduction of this cut- 
off wavenumber removes the divergence in the integral, 
but it also changes the above scaling argument. Since the 
dominant contribution to Eq. (C3) still comes from small 
q, we need to identify the scaling of the mode-coupling 
coefficients with A: for (7 <C k, M\k,-\k.\q — A™ M^.-k.q- 
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Cascade 

2D enstrophy 

2D energy 

3D energy 

3D helicity 



nDiA 



TABLE I. Scaling exponents for various cascades in two 
dimensions (2D) and three dimensions (3D), using either the 
streamfunction = $ or velocity ^1) = u normalization. 



Since the lower wavenumber limit is now fixed, no self- 
similar scaling in q can be made; the scaling with k for 
small q then leads to \^'^+^^' = l. Hence for the DIA 
equations the actual scalings of the response function, 
correlation function, and energy spectrum are given by 



^DIA 



-m 



nDiA 



/3. 



DIA 



7 ■ 



m 



1 — m 



(C8a) 



(C8b) 



(C8c) 



In Table | we compare the scalings in Eqs . (C7) with 
the anomalous DI A sc alings given by Eq. ( |C8|) . The scal- 
ings given by Eq. (C7) are consistent with Kolmogorov's 
dimensional analysis. We emphasize that these scal- 
ings would have also been obtained for the DIA equa- 
tions (they too are dimensionally consistent) had the 
wavenumber integral in Eq. (C5) converged. 
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